xi_func=@(rho)(1-rq_ratio)*max([1+rho,1/(1+pi_m)])+rq_ratio/(1+pi_m)-cost;
rho_bar = fzero(@(x)xi_func(x)-1/beta,[1e-10,2]);
rho_bar = min(rho_bar,1/beta-1+cost);
rho_grid=linspace(1/(1+pi_m)-1,rho_bar,100);

func_var.beta=beta;
func_var.alpha1=alpha1;
func_var.alpha2=alpha2;
func_var.alpha3=alpha3;
func_var.lambda_func=lambda_func;
func_var.dlambda_func=dlambda_func;
func_var.rq_ratio=rq_ratio;
func_var.cost=cost;
func_var.df_inv = df_inv;
func_var.df=df;
func_var.p_bar=p_grid(end);
func_var.rho_bar = rho_bar;
func_var.Nb=Nb;
func_var.i_m=i_m;